Técnicas de clasificación y evaluación de procesos en sistemas forestales

Jessica Bernal Borrego

Contenido:

  1. Introducción
  2. Creación muestreo
  3. Clasficación supervisada
  4. Clasificador KNN
  5. Clasificador ANN
  6. Clasificador SVM
  7. Clasificador RF
  8. comparativa resultados entre clasificadores
  9. Clasificador multi-temporal

###Introducción

Los objetivos de este tutorial son: - Clasificar una extracto de una escena Sentinel 2 usando diferentes clasificadores de machine learning. - Comparar los resultados obtenidos empleando distintos clasificadores

Los paquete que emplearemos son:

library(sp)
library(rgdal)
library(raster)
library(reshape)
library(grid)
library(gridExtra)
library(RStoolbox)
library(caret)
library(rasterVis)
library(corrplot)
library(doParallel)
library(NeuralNetTools)
library(tidyr)
library(stringr)
library(e1071)
library(sf)
library (mapview) 

El siguiente paso consistirá en definir el directorio de trabajo donde se localizarán nuestras imágenes. Los datos se corresponden con una escena Sentinel 2 (fichero denominado sentinel_o_bxx.tif, siendo xx el número de la banda)

setwd("C:/Geoforest/Tec_Clasif")
dir_in<-"./Material_practicas/Sentinel/O"

Como en cuadernos anteriores creamos una lista con los nombres de los archivos alojados en el directorio de trabajo, creando posteriormente un rasterstack

rasList <- list.files(dir_in,pattern="tif",
                      full.names = TRUE)

sentinel_o <- stack(rasList)

Una vez generado, vamos a proceder a comprobar los atributos de la escena

sentinel_o
## class      : RasterStack 
## dimensions : 2052, 2057, 4220964, 10  (nrow, ncol, ncell, nlayers)
## resolution : 20, 20  (x, y)
## extent     : 300920, 342060, 4042100, 4083140  (xmin, xmax, ymin, ymax)
## crs        : +proj=utm +zone=30 +datum=WGS84 +units=m +no_defs 
## names      : sentinel_o_b01, sentinel_o_b02, sentinel_o_b03, sentinel_o_b04, sentinel_o_b05, sentinel_o_b06, sentinel_o_b07, sentinel_o_b08, sentinel_o_b09, sentinel_o_b10 
## min values :              0,              0,              0,              0,              0,              0,              0,              0,              0,              0 
## max values :          65535,          65535,          65535,          65535,          65535,          65535,          65535,          65535,          65535,          65535

A continuación se van a generar el gráfico de densidad y el histograma para una de las bandas. En primer lugar lo calcularemos para una de las bandas (imagen).

gdensidad=ggplot(sentinel_o,aes(sentinel_o_b01))+geom_density()
ghistograma=ggplot(sentinel_o,aes(sentinel_o_b01))+geom_histogram()

print (gdensidad)

print (ghistograma)
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Si las queremos hacer de todas las bandas es posible generar un gráfico por cada una de ellas y posteriormente componerlos en una sola figura

gdens1=ggplot(sentinel_o,aes(sentinel_o_b01))+geom_density()
gdens2=ggplot(sentinel_o,aes(sentinel_o_b02))+geom_density()
gdens3=ggplot(sentinel_o,aes(sentinel_o_b03))+geom_density()
gdens4=ggplot(sentinel_o,aes(sentinel_o_b04))+geom_density()
gdens5=ggplot(sentinel_o,aes(sentinel_o_b05))+geom_density()
gdens6=ggplot(sentinel_o,aes(sentinel_o_b06))+geom_density()
gdens7=ggplot(sentinel_o,aes(sentinel_o_b07))+geom_density()
gdens8=ggplot(sentinel_o,aes(sentinel_o_b08))+geom_density()
gdens9=ggplot(sentinel_o,aes(sentinel_o_b09))+geom_density()
gdens10=ggplot(sentinel_o,aes(sentinel_o_b10))+geom_density()

grid.arrange(gdens1,gdens2,gdens3,gdens4,gdens5,gdens6,gdens7,gdens8,gdens9,gdens10,ncol=4,nrow=4)

Creación muestreo

Al ser una clasificación supervisada necesitaremos aportar al clasificador la información necesaria para realizar las fases de entrenamiento y validación. A partir del MFE facilitado, con geometría poligonal, tenemos dos opciones a la hora de continuar. Para ello, es necesario saber que en este entorno de trabajo la información suministrada al clasificador deberá ser de tipo poligonal. Por ello, las opciones son:

  • Opción 1: Preparación de los datos desde un software externo, por ejemplo QGIS, y luego leerlo en R.

  • Opción 2: Preparación de los datos desde R.

En el caso de optar por trabajar con una herramienta externa los datos de tipo puntual pueden ser almacendados en formato shapefile y luego ser leidos mediante la función readOGR.

A modo de ejemplo, la llamada a la función sería:

#train_data<-readOGR('./entrenamiento/muestreo_500.shp')

Debiendo comprobar posteriormente la estructura del fichero leido.

#str(train_data)

Tambien es posible hacer la misma operación mediante la función st_read.

En este cuaderno se optará por trabajar con la segunda opción, es decir, preparando el muestreo desde R.

Para ello, en primer lugar se procederá a leer el archivo shapefile del MFE.

#Generamos una semilla para garantizar la repetitividad de los resultados
set.seed(123)
MFE=st_read('./Material_practicas/MFE/MFE.shp')
## Reading layer `MFE' from data source 
##   `C:\Geoforest\Tec_Clasif\Material_practicas\MFE\MFE.shp' using driver `ESRI Shapefile'
## Simple feature collection with 22 features and 5 fields
## Geometry type: MULTIPOLYGON
## Dimension:     XY
## Bounding box:  xmin: 310927.3 ymin: 4052098 xmax: 332056.3 ymax: 4073143
## Projected CRS: ETRS89 / UTM zone 30N
mapview (MFE,zcol='leyenda')

Con objeto de obtener una muestra balanceada se va a determinar el total de la superficie ocupada por cada clase de la leyenda, repartiendo el tamaño de la muestra proporcional a la superficie ocupada.

clases = unique(MFE$leyenda)

area_total = sum(st_area(MFE))
area_clases=0
for (i in 1:length(clases)) {
  geom_clase=MFE[which(MFE$leyenda == clases[i],arr.ind=FALSE),]
  area_clases[i]=sum(st_area(geom_clase))
}

Se va a fijar un tamaño total de muestreo igual a 500 puntos de tal forma que en cada clase habrá el siguiente número de muestras.

num_muestras= as.integer(500*area_clases/area_total)
print(num_muestras)
##  [1]  28   0  26   6   0   0  15   0  25  36 122  34 203

En este punto, analizar el número de cada clase. ¿Hay muestras en todas las clases?, ¿Que significa que haya clases con un número elevado y otras que sea muy reducido o directamente cero?

Aun sabiendo que no es correcto, en lugar de establecer el muestreo atendiendo al criterio anterior se va a seleccionar un número fijo de muestras para todas las clases. De esta manera, el trabajo a entregar por el alumno será determinar una leyenda adecuada a la variabilidad espacial y espectral de las clases presentes en la escena.

Por ello, en primer lugar se va a proceder a realizar un muestro sobre el MFE de tipo aleatorio, extrayendo la información temática a partir de la función st_join.

puntos.ref <- st_sample(MFE, c(50,50,50,50,50,50,50,50,50,50,50,50,50), type='random',exact=TRUE) #Generamos una lista de puntos de forma aleatoria
puntos.ref<-st_sf(puntos.ref) #Convertimos la lista en un spatial feature

puntos.ref<-st_join(puntos.ref,MFE) #Cruzamos los datos
puntos.ref_backup <- puntos.ref
mapview(puntos.ref, zcol='leyenda')

Y su representación sobre el MFE:

mapview(MFE, zcol='leyenda')+mapview(puntos.ref,zcol='leyenda')

A continuación se va a proceder a obtener la firma espectral de cada una de las clases. en primer lugar será necesario extraer los valores de reflectancia para cada punto en cada una de las bandas mediante el comando extract. Como estos datos los representaremos mediante la librería ggplot los convertiremos a un tipo de dato dataframe. Además, antes de la representación se determinarán los valores medios de reflectancia por clase y banda.

Así, en primer lugar generamos el dataframe.

puntos.ref=as_Spatial(puntos.ref)
puntos.ref@data$leyenda=as.factor(puntos.ref@data$leyenda)

reflectancia<- as.data.frame(raster::extract(sentinel_o,puntos.ref))

Comprobando los valores extraidos.

head(reflectancia)
##   sentinel_o_b01 sentinel_o_b02 sentinel_o_b03 sentinel_o_b04 sentinel_o_b05
## 1            935            701            434            684           1405
## 2           1069            871            703            985           2058
## 3            869            562            356            456            800
## 4            875            583            360            454            843
## 5            922            681            469            694           1312
## 6            929            696            471            702           1611
##   sentinel_o_b06 sentinel_o_b07 sentinel_o_b08 sentinel_o_b09 sentinel_o_b10
## 1           1705           1545           1754            689            340
## 2           2270           2263           2531           1983           1174
## 3            962            787            883            497            258
## 4            853            901           1032            755            333
## 5           1720           1605           1697           1248            698
## 6           1998           1923           2179           1293            633

Calculamos el valor medio de reflectancia de cada clase para cada banda. Para ello, usaremos la función aggregate() para unir los puntos de entrenamiento por clase.

mean_reflectancia <-aggregate(reflectancia,list(puntos.ref$leyenda),mean,na.rm = TRUE)

Comprobamos los valores medios obtenidos

head(mean_reflectancia)
##                 Group.1 sentinel_o_b01 sentinel_o_b02 sentinel_o_b03
## 1                    --       1088.940         954.56         782.28
## 2        Bosques Mixtos       1040.033         842.18         704.48
## 3            Cadufolios       1348.460        1171.52        1135.10
## 4 Castaños/Caducifolios       1105.320         954.24        1059.78
## 5             Eucalipto        926.980         672.34         475.60
## 6              Matorral       1460.080        1292.59        1219.90
##   sentinel_o_b04 sentinel_o_b05 sentinel_o_b06 sentinel_o_b07 sentinel_o_b08
## 1        1126.10       2144.880       2429.600        2449.42        2619.32
## 2         910.86       1535.053       1722.407        1719.52        1846.22
## 3        1327.40       1725.660       1868.640        1904.24        2057.20
## 4        1233.26       1441.640       1580.820        1616.80        1753.70
## 5         629.52       1082.500       1208.240        1133.80        1279.34
## 6        1449.83       2027.840       2184.040        2173.93        2315.40
##   sentinel_o_b09 sentinel_o_b10
## 1       1576.540         945.54
## 2       1255.907         784.16
## 3       1897.600        1326.12
## 4       2244.360        1609.96
## 5        784.000         439.12
## 6       1613.840        1137.70

Por la forma en la que estan almacenados los datos en el dataframe (cada banda se almacena en una columna) es necesario modificarlo para que esten todos los datos de reflectancias registrados en una columna, creando una nueva columna donde se registre la banda de donde proceden, de tal forma que la información aparecerá ordenada por filas.

mean_reflectance2 <- gather(mean_reflectancia, key="banda", value="reflectance",sentinel_o_b01:sentinel_o_b10)

Si analizamos el contenido del dataframe no se dispone de un campo numerico que permita ordenar las bandas a la hora de pintarlas. Por ello se va a crear un nuevo campo de tipo numérico con el número de la banda.

mean_reflectance2$banda_num=(as.numeric(str_replace(mean_reflectance2$banda,"sentinel_o_b","")))

Finalmente, mediante ggplot se pintará un gráfico de tipo geom_line para ver la firma espectral de cada una de clases.

Como se puede ver, muchas de las clases presentan un comportamiento muy similar y por tanto la calidad temática de los resultados de la clasificación a priori pueden ser bajos.

ggplot(mean_reflectance2,aes(x=banda_num,y=reflectance))+
  geom_line(aes(colour = Group.1))+theme_bw()

Clasficación supervisada

En primer lugar veremos el uso de alguno de los operadores clásicos empleados en Teledetección para clasificar imágenes, en este caso clasificador por máxima probabilidad. Para ello se empleará la función superClass dentro del paquete RStoolbox. De entre los parámetros a incluir en la función destacar que:

  • trainData contiene el muestreo a emplear en la clasificación.

  • trainPartition: contiene un valor indicando el tamaño de la muestra destinada al entrenamiento.

  • model indica el tipo de clasficación que se desea realizar, por defecto la función aplica un random forest (rf) pero en este caso se realizará una clasificación por máxima probabilidad (mlc).

#puntos.ref<- as_Spatial(puntos.ref)
puntos.ref@data$leyenda=as.factor(puntos.ref@data$leyenda)

Max.prob<- superClass(sentinel_o, 
                      trainData = puntos.ref, 
                      trainPartition =0.5,
                      responseCol = "leyenda",
                      model = "mlc") #máxima probabilidad

El resultado de la clasificación se muestra recogido en una variable de tipo lista. En el quinto elemento de la lista se recoge el resultado cartográfico mediante un rasterLayer, pudiendo ser representado por ejemplo mediante la función plot.

leyenda_colores <- viridis::viridis(13)
plot(Max.prob$map,
     col=leyenda_colores,
     legend = FALSE)
legend("topright",cex=0.65, y.intersp = 0.55,x.intersp = 0.5,
        legend = levels(as.factor(puntos.ref$leyenda)),
        fill = leyenda_colores ,title = "",
        inset=c(0,0))

Además del producto cartográfico es necesario realizar una evaluación de la calidad temática. Para ello mediante en el segundo elemento de la lista, denominado modelFit se encuentra la matriz de confusión resultante así como los valores de exactitud global y kapp obtenidos en el entrenamiento. Por otro lado en el elemento results aparecen recogidos estos elementos de calidad global y su desviación.

Max.prob$modelFit
## [[1]]
##   TrainAccuracy TrainKappa method
## 1      0.354396  0.2797632 custom
## 
## [[2]]
## Cross-Validated (5 fold) Confusion Matrix 
## 
## (entries are average cell counts across resamples)
##  
##                        Reference
## Prediction                -- Bosques Mixtos Cadufolios Castaños/Caducifolios
##   --                     1.8            0.2        0.0                   0.0
##   Bosques Mixtos         0.2            3.8        1.0                   0.0
##   Cadufolios             0.0            0.0        1.0                   0.4
##   Castaños/Caducifolios  0.0            0.0        0.0                   3.8
##   Eucalipto              0.2            1.4        0.2                   0.0
##   Matorral               0.2            0.2        0.4                   0.0
##   Pinos                  1.6            3.0        0.0                   0.0
##   Pinos/Pinsapos         0.2            0.8        0.0                   0.0
##   Pinsapos               0.0            0.4        0.0                   0.0
##   Quercineas             0.0            2.6        0.4                   0.0
##   Ribera                 0.0            0.6        0.0                   0.0
##   Suelos?                0.2            0.8        1.0                   0.2
##   Veg. Dispersa          0.2            0.8        1.0                   0.4
##                        Reference
## Prediction              Eucalipto Matorral Pinos Pinos/Pinsapos Pinsapos
##   --                          0.0      0.0   1.2            0.6      0.0
##   Bosques Mixtos              0.2      1.6   0.6            0.6      0.2
##   Cadufolios                  0.0      0.4   0.0            0.0      0.2
##   Castaños/Caducifolios       0.0      0.0   0.0            0.0      0.0
##   Eucalipto                   1.8      0.6   0.6            0.4      0.2
##   Matorral                    0.2      3.6   0.0            0.0      0.2
##   Pinos                       0.6      0.8  10.8            1.0      1.0
##   Pinos/Pinsapos              1.2      0.8   2.0            1.0      0.8
##   Pinsapos                    0.2      0.4   2.2            0.4      2.0
##   Quercineas                  0.6      0.4   1.4            1.0      0.0
##   Ribera                      0.0      0.0   0.2            0.0      0.2
##   Suelos?                     0.0      0.0   0.6            0.0      0.2
##   Veg. Dispersa               0.2      1.4   0.4            0.0      0.0
##                        Reference
## Prediction              Quercineas Ribera Suelos? Veg. Dispersa
##   --                           0.2    0.0     0.2           0.2
##   Bosques Mixtos               0.8    0.2     1.4           1.8
##   Cadufolios                   0.2    0.0     0.2           1.6
##   Castaños/Caducifolios        0.0    0.0     0.0           0.0
##   Eucalipto                    0.8    0.2     0.0           0.4
##   Matorral                     0.0    0.2     0.0           1.6
##   Pinos                        0.8    1.2     1.4           1.6
##   Pinos/Pinsapos               0.2    0.2     0.0           0.6
##   Pinsapos                     0.0    0.6     0.4           0.2
##   Quercineas                   6.6    0.8     0.0           4.2
##   Ribera                       0.4    0.2     0.0           0.0
##   Suelos?                      0.0    0.0     0.4           1.4
##   Veg. Dispersa                0.0    0.2     1.0           1.4
##                             
##  Accuracy (average) : 0.3544
Max.prob$model$results
##   parameter Accuracy     Kappa AccuracySD    KappaSD
## 1      none 0.354396 0.2797632 0.02372322 0.02503495

Creación de un dataframe con los puntos del entrenamiento etiquetados y con sus valores de reflectancia

En el paso anterior se generó un dataframe con los valores medios de reflectancia para cada clase. Ahora se va a generar un dataframe que contendrá para cada punto su etiqueta y los valores de reflectancia de todas y cada una de las bandas. Advertir que el objeto apunta a toda esa información, de forma que *@* es un operador especial que permite acceder a un objeto dentro de otro objeto. Nota: Se va a crear una variable denominada train_data igual a puntos.ref por si el trascurso de la práctica se comete un error, pudiendo tener así un backup de los datos hasta este punto.

train_data = as_Spatial(puntos.ref_backup)
train_data@data$leyenda=as.factor(puntos.ref@data$leyenda)
train_data@data=data.frame(train_data@data,reflectancia[match(rownames(train_data@data),rownames(reflectancia)),])

Veamos el resultado obtenido.

str(train_data)
## Formal class 'SpatialPointsDataFrame' [package "sp"] with 5 slots
##   ..@ data       :'data.frame':  1100 obs. of  15 variables:
##   .. ..$ OBJECTID      : num [1:1100] 1 1 1 1 1 1 1 1 1 1 ...
##   .. ..$ NOM_FORARB    : chr [1:1100] "Alcornocales" "Alcornocales" "Alcornocales" "Alcornocales" ...
##   .. ..$ Shape_Leng    : num [1:1100] 34046 34046 34046 34046 34046 ...
##   .. ..$ Shape_Area    : num [1:1100] 4419619 4419619 4419619 4419619 4419619 ...
##   .. ..$ leyenda       : Factor w/ 13 levels "--","Bosques Mixtos",..: 10 10 10 10 10 10 10 10 10 10 ...
##   .. ..$ sentinel_o_b01: num [1:1100] 935 1069 869 875 922 ...
##   .. ..$ sentinel_o_b02: num [1:1100] 701 871 562 583 681 696 609 753 500 654 ...
##   .. ..$ sentinel_o_b03: num [1:1100] 434 703 356 360 469 471 377 542 282 391 ...
##   .. ..$ sentinel_o_b04: num [1:1100] 684 985 456 454 694 702 477 791 254 562 ...
##   .. ..$ sentinel_o_b05: num [1:1100] 1405 2058 800 843 1312 ...
##   .. ..$ sentinel_o_b06: num [1:1100] 1705 2270 962 853 1720 ...
##   .. ..$ sentinel_o_b07: num [1:1100] 1545 2263 787 901 1605 ...
##   .. ..$ sentinel_o_b08: num [1:1100] 1754 2531 883 1032 1697 ...
##   .. ..$ sentinel_o_b09: num [1:1100] 689 1983 497 755 1248 ...
##   .. ..$ sentinel_o_b10: num [1:1100] 340 1174 258 333 698 ...
##   ..@ coords.nrs : num(0) 
##   ..@ coords     : num [1:1100, 1:2] 326741 329242 328773 331246 329865 ...
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : NULL
##   .. .. ..$ : chr [1:2] "lon" "lat"
##   ..@ bbox       : num [1:2, 1:2] 312560 4052227 331754 4072825
##   .. ..- attr(*, "dimnames")=List of 2
##   .. .. ..$ : chr [1:2] "lon" "lat"
##   .. .. ..$ : chr [1:2] "min" "max"
##   ..@ proj4string:Formal class 'CRS' [package "sp"] with 1 slot
##   .. .. ..@ projargs: chr "+proj=utm +zone=30 +ellps=GRS80 +towgs84=0,0,0,0,0,0,0 +units=m +no_defs"
##   .. .. ..$ comment: chr "PROJCRS[\"ETRS89 / UTM zone 30N\",\n    BASEGEOGCRS[\"ETRS89\",\n        DATUM[\"European Terrestrial Reference"| __truncated__

Podemos observar como nos indica que tenemos 11 variables, una correspondiente a la clase y diez a las bandas espectrales, siendo estas últimas nuestras variables predictoras de la variable respuesta consistente en las clases. Ahora, podriamos ver un resumen estadístico de la variable.

summary(train_data@data)
##     OBJECTID     NOM_FORARB          Shape_Leng       Shape_Area      
##  Min.   : 1.0   Length:1100        Min.   :  1805   Min.   :   52512  
##  1st Qu.: 6.0   Class :character   1st Qu.:  3900   1st Qu.:   84213  
##  Median :11.5   Mode  :character   Median : 26445   Median : 2351264  
##  Mean   :11.5                      Mean   : 64448   Mean   : 9131237  
##  3rd Qu.:17.0                      3rd Qu.: 80456   3rd Qu.:10093586  
##  Max.   :22.0                      Max.   :455519   Max.   :81665992  
##                                                                       
##            leyenda    sentinel_o_b01 sentinel_o_b02   sentinel_o_b03  
##  Pinos         :200   Min.   : 761   Min.   : 457.0   Min.   : 249.0  
##  Bosques Mixtos:150   1st Qu.: 888   1st Qu.: 632.8   1st Qu.: 408.0  
##  Veg. Dispersa :150   Median : 980   Median : 787.5   Median : 587.5  
##  Matorral      :100   Mean   :1103   Mean   : 905.1   Mean   : 775.2  
##  Quercineas    :100   3rd Qu.:1126   3rd Qu.: 995.2   3rd Qu.: 963.2  
##  --            : 50   Max.   :5384   Max.   :5290.0   Max.   :5925.0  
##  (Other)       :350                                                   
##  sentinel_o_b04   sentinel_o_b05 sentinel_o_b06 sentinel_o_b07 sentinel_o_b08
##  Min.   : 226.0   Min.   : 247   Min.   : 233   Min.   : 212   Min.   : 187  
##  1st Qu.: 561.5   1st Qu.:1030   1st Qu.:1158   1st Qu.:1146   1st Qu.:1234  
##  Median : 849.0   Median :1566   Median :1770   Median :1766   Median :1918  
##  Mean   : 979.8   Mean   :1561   Mean   :1737   Mean   :1732   Mean   :1869  
##  3rd Qu.:1241.5   3rd Qu.:2002   3rd Qu.:2242   3rd Qu.:2275   3rd Qu.:2448  
##  Max.   :6146.0   Max.   :6331   Max.   :6394   Max.   :6354   Max.   :6308  
##                                                                              
##  sentinel_o_b09   sentinel_o_b10  
##  Min.   :  53.0   Min.   :  28.0  
##  1st Qu.: 600.5   1st Qu.: 309.2  
##  Median :1121.0   Median : 639.0  
##  Mean   :1300.0   Mean   : 834.0  
##  3rd Qu.:1829.5   3rd Qu.:1186.0  
##  Max.   :5973.0   Max.   :4484.0  
## 

En este caso se observa como no hay datos ausentes (NA). En caso de que aparecieran es importante eliminarlos, empleando para ello la función na.omit(). Una vez aplicada podemos emplear la función complete.cases() para comprobar que se han borrado.

train_data@data= na.omit(train_data@data)
complete.cases(train_data@data)
##    [1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##   [15] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##   [29] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##   [43] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##   [57] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##   [71] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##   [85] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##   [99] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [113] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [127] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [141] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [155] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [169] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [183] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [197] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [211] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [225] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [239] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [253] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [267] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [281] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [295] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [309] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [323] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [337] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [351] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [365] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [379] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [393] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [407] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [421] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [435] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [449] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [463] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [477] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [491] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [505] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [519] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [533] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [547] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [561] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [575] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [589] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [603] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [617] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [631] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [645] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [659] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [673] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [687] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [701] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [715] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [729] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [743] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [757] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [771] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [785] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [799] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [813] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [827] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [841] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [855] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [869] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [883] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [897] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [911] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [925] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [939] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [953] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [967] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [981] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
##  [995] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1009] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1023] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1037] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1051] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1065] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1079] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
## [1093] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE

Preparación del set de entrenamiento y testeo

Es recomendable separar de forma aleatoria el set de entrenamiento inicial en 3 grupos: entrenamiento, validación y testeo. En este caso solo lo vamos a separar en los dos primeros. Para ello, en primer lugar es necesario establecer un valor predefinido de semilla empleando set.seed().

hre_seed<- 123
set.seed(hre_seed)

Ahora, dividiremos nuestro set de entrenamiento inicial en entrenamiento y testeo usando para ello la función createDataPartition() del paquete caret. Por ejemplo, vamos a establecer que el 80% de los datos iniciales pasen a ser de entrenamiento y el 20% de test. Recordar que el entrenamiento nos permitirá optimizar los parámetros iniciales del modelo mientras que los de testeo nos permitirán evaluar la calidad del mismo.

Nota: Se ha establecido como variable $leyenda pues es esta la que que contiene la etiqueta de las clases. Por otro lado el parámetro p contiene el porcentaje a emplear en la separación de la muestra. Finalmente, el parámetro list indica si devuelve una lista o una matriz, en nuestro caso indicaremos FALSE de forma que devuelve una matriz.

Así, el resultado de la variable inTraining como podemos comprobar es una lista de valores núméricos indicando el índice de los elementos empleados para entrenamiento.

inTraining <- createDataPartition(train_data@data$leyenda, p=0.80,list=FALSE)

training <-train_data@data[inTraining,]
training=training[,-(1:4)] #Borramos columnas no necesarias

testing <- train_data@data[-inTraining,]
testing=testing[,-(1:4)] #Borramos columnas no necesarias

Resumen estadistico de los set de entreneamiento y testeo

Antes de comenzar el entrenamiento del clasificador de machine learning previo a la clasificación es necesario realizar un chequeo de los set de datos pues puede ser que las imágenes presenten problemas o que hayamos cometido errores en la identificación. Así, en primer lugar vamos a obtener un resumen estadístico de ambos set.

summary(training)
##            leyenda    sentinel_o_b01   sentinel_o_b02   sentinel_o_b03  
##  Pinos         :160   Min.   : 761.0   Min.   : 457.0   Min.   : 249.0  
##  Bosques Mixtos:120   1st Qu.: 884.0   1st Qu.: 624.0   1st Qu.: 402.0  
##  Veg. Dispersa :120   Median : 983.5   Median : 792.0   Median : 593.5  
##  Matorral      : 80   Mean   :1100.7   Mean   : 903.2   Mean   : 775.1  
##  Quercineas    : 80   3rd Qu.:1126.0   3rd Qu.: 995.8   3rd Qu.: 973.8  
##  --            : 40   Max.   :5342.0   Max.   :5290.0   Max.   :5925.0  
##  (Other)       :280                                                     
##  sentinel_o_b04   sentinel_o_b05 sentinel_o_b06 sentinel_o_b07 sentinel_o_b08
##  Min.   : 226.0   Min.   : 247   Min.   : 233   Min.   : 212   Min.   : 187  
##  1st Qu.: 543.0   1st Qu.:1017   1st Qu.:1137   1st Qu.:1114   1st Qu.:1204  
##  Median : 862.0   Median :1560   Median :1771   Median :1774   Median :1924  
##  Mean   : 978.8   Mean   :1555   Mean   :1729   Mean   :1723   Mean   :1862  
##  3rd Qu.:1243.0   3rd Qu.:2010   3rd Qu.:2250   3rd Qu.:2278   3rd Qu.:2443  
##  Max.   :6146.0   Max.   :6331   Max.   :6394   Max.   :6354   Max.   :6308  
##                                                                              
##  sentinel_o_b09 sentinel_o_b10  
##  Min.   :  53   Min.   :  28.0  
##  1st Qu.: 580   1st Qu.: 301.0  
##  Median :1138   Median : 646.5  
##  Mean   :1315   Mean   : 845.8  
##  3rd Qu.:1848   3rd Qu.:1196.5  
##  Max.   :5973   Max.   :4484.0  
## 
summary(testing)
##            leyenda   sentinel_o_b01   sentinel_o_b02   sentinel_o_b03  
##  Pinos         :40   Min.   : 783.0   Min.   : 496.0   Min.   : 282.0  
##  Bosques Mixtos:30   1st Qu.: 899.5   1st Qu.: 662.0   1st Qu.: 434.0  
##  Veg. Dispersa :30   Median : 972.0   Median : 772.5   Median : 568.0  
##  Matorral      :20   Mean   :1110.8   Mean   : 912.7   Mean   : 775.9  
##  Quercineas    :20   3rd Qu.:1123.8   3rd Qu.: 995.2   3rd Qu.: 917.5  
##  --            :10   Max.   :5384.0   Max.   :5289.0   Max.   :5881.0  
##  (Other)       :70                                                     
##  sentinel_o_b04   sentinel_o_b05 sentinel_o_b06 sentinel_o_b07 sentinel_o_b08
##  Min.   : 264.0   Min.   : 303   Min.   : 267   Min.   : 253   Min.   : 213  
##  1st Qu.: 612.5   1st Qu.:1154   1st Qu.:1264   1st Qu.:1254   1st Qu.:1335  
##  Median : 807.0   Median :1576   Median :1751   Median :1749   Median :1904  
##  Mean   : 983.7   Mean   :1588   Mean   :1766   Mean   :1767   Mean   :1900  
##  3rd Qu.:1219.2   3rd Qu.:1958   3rd Qu.:2195   3rd Qu.:2262   3rd Qu.:2457  
##  Max.   :5873.0   Max.   :5905   Max.   :5943   Max.   :6284   Max.   :6112  
##                                                                              
##  sentinel_o_b09   sentinel_o_b10  
##  Min.   :  87.0   Min.   :  47.0  
##  1st Qu.: 645.8   1st Qu.: 320.0  
##  Median :1094.0   Median : 596.5  
##  Mean   :1241.2   Mean   : 787.0  
##  3rd Qu.:1681.2   3rd Qu.:1072.5  
##  Max.   :4716.0   Max.   :3370.0  
## 

Posteriormente vamos a generar un grafico de densidades para cada clase / banda que permita representar la distribución de los datos. Esto va a permitirnos evaluar si hay una adecuada separabilidad entre las clases. Además permite determinar si el efecto cizalla en la distribución es acusado o no. Nota: Se han seleccionado los índices 2 al 11 pues en el caso de este ejemplo contienen los datos de reflectancia para cada punto en todas las bandas empleadas.

featurePlot(x=training[,2:11],
            y=training$leyenda,
            plot="density",
            labels=c("Reflectancia","Distribucion densidades"),
            layout=c(2,2))

Por otro lado podemos calcular la cizalla mediante la función skewness() de la librería e1071. Nos apartorá información si al distribución es simétrica o no. Por lo general, una distribución es simétrica cuando el valor de skewness es 0 o próximo a 0.

skwenessvalues <- apply(training[,2:11],2,skewness)
skwenessvalues
## sentinel_o_b01 sentinel_o_b02 sentinel_o_b03 sentinel_o_b04 sentinel_o_b05 
##      4.6181023      3.6657509      3.0889411      2.4448703      0.8939717 
## sentinel_o_b06 sentinel_o_b07 sentinel_o_b08 sentinel_o_b09 sentinel_o_b10 
##      0.5272516      0.4532052      0.3274439      1.1063450      1.4029482

Por otro lado, si se detecta alguna distribución bimodal puede ser indicativo de una posible presencia de errores groseros en el muestreo. De forma complementaria podemos representar los datos mediante cajas de bigotes con objeto de ver esta presencia.

featurePlot(x=training[,2:11],
            y=training$leyenda,
            plot="box",
            layout=c(2,2))

La posible presencia de errores groseros podría deberse por una parte a fallos humanos o tambien a la propia variabilidad del territorio y el comportamiento de las clases. Supongamos por ejemplo que contamos con las clases “uso agricola” y “suelo desnudo”, es posible que estas dos clases presenten un comportamiento similar, siendo adecuado analizar la correlación entre bandas.

A modo de ejemplo se presentan graficos de correlación entre dos clases y 6 bandas espectrales, viendo una clara correlación entre bandas.

band1_2 <-ggplot(data=training,aes(sentinel_o_b01,sentinel_o_b02))+
                   geom_point(aes(shape=leyenda,colour=leyenda))

band1_3 <-ggplot(data=training,aes(sentinel_o_b01,sentinel_o_b03))+
                   geom_point(aes(shape=leyenda,colour=leyenda))

band1_4 <-ggplot(data=training,aes(sentinel_o_b01,sentinel_o_b04))+
                   geom_point(aes(shape=leyenda,colour=leyenda))

band1_5 <-ggplot(data=training,aes(sentinel_o_b01,sentinel_o_b05))+
                   geom_point(aes(shape=leyenda,colour=leyenda))

grid.arrange(band1_2,band1_3,band1_4,band1_5)
## Warning: The shape palette can deal with a maximum of 6 discrete values because
## more than 6 becomes difficult to discriminate; you have 13. Consider
## specifying shapes manually if you must have them.
## Warning: Removed 520 rows containing missing values (geom_point).
## Warning: The shape palette can deal with a maximum of 6 discrete values because
## more than 6 becomes difficult to discriminate; you have 13. Consider
## specifying shapes manually if you must have them.
## Warning: Removed 520 rows containing missing values (geom_point).
## Warning: The shape palette can deal with a maximum of 6 discrete values because
## more than 6 becomes difficult to discriminate; you have 13. Consider
## specifying shapes manually if you must have them.
## Warning: Removed 520 rows containing missing values (geom_point).
## Warning: The shape palette can deal with a maximum of 6 discrete values because
## more than 6 becomes difficult to discriminate; you have 13. Consider
## specifying shapes manually if you must have them.
## Warning: Removed 520 rows containing missing values (geom_point).

Numéricamente, mediante la función cor() calculamos la correlacion entre las bandas espectrales de la escena. El resultado puede ser “complejo” de analizar, siendo mejor una representación gráfica.

bandcorrelaciones = cor(training[,2:11])

bandcorrelaciones
##                sentinel_o_b01 sentinel_o_b02 sentinel_o_b03 sentinel_o_b04
## sentinel_o_b01      1.0000000      0.9859977      0.9640630      0.9285340
## sentinel_o_b02      0.9859977      1.0000000      0.9882650      0.9734298
## sentinel_o_b03      0.9640630      0.9882650      1.0000000      0.9811764
## sentinel_o_b04      0.9285340      0.9734298      0.9811764      1.0000000
## sentinel_o_b05      0.7511718      0.8281299      0.8091457      0.8948271
## sentinel_o_b06      0.6895358      0.7743948      0.7554891      0.8514766
## sentinel_o_b07      0.6781686      0.7679473      0.7539346      0.8491108
## sentinel_o_b08      0.6479488      0.7417461      0.7300217      0.8323305
## sentinel_o_b09      0.4105274      0.5346405      0.6079181      0.6789671
## sentinel_o_b10      0.4526035      0.5684750      0.6488496      0.6984570
##                sentinel_o_b05 sentinel_o_b06 sentinel_o_b07 sentinel_o_b08
## sentinel_o_b01      0.7511718      0.6895358      0.6781686      0.6479488
## sentinel_o_b02      0.8281299      0.7743948      0.7679473      0.7417461
## sentinel_o_b03      0.8091457      0.7554891      0.7539346      0.7300217
## sentinel_o_b04      0.8948271      0.8514766      0.8491108      0.8323305
## sentinel_o_b05      1.0000000      0.9912350      0.9877215      0.9837551
## sentinel_o_b06      0.9912350      1.0000000      0.9936050      0.9932052
## sentinel_o_b07      0.9877215      0.9936050      1.0000000      0.9938419
## sentinel_o_b08      0.9837551      0.9932052      0.9938419      1.0000000
## sentinel_o_b09      0.6890304      0.6935422      0.7193895      0.7380472
## sentinel_o_b10      0.6491981      0.6426895      0.6703054      0.6839346
##                sentinel_o_b09 sentinel_o_b10
## sentinel_o_b01      0.4105274      0.4526035
## sentinel_o_b02      0.5346405      0.5684750
## sentinel_o_b03      0.6079181      0.6488496
## sentinel_o_b04      0.6789671      0.6984570
## sentinel_o_b05      0.6890304      0.6491981
## sentinel_o_b06      0.6935422      0.6426895
## sentinel_o_b07      0.7193895      0.6703054
## sentinel_o_b08      0.7380472      0.6839346
## sentinel_o_b09      1.0000000      0.9862417
## sentinel_o_b10      0.9862417      1.0000000

Esta sería la representación gráfica de la matriz de correlación

corrplot(bandcorrelaciones,method="number")

corrplot(bandcorrelaciones,method="number",type = "upper")

corrplot(bandcorrelaciones,method="color",type="lower")

Definición de los parámetros del modelo para entrenamiento

Este paso es uno de los mas importantes, pued de la correcta configuración de los parámetros dependerán nuestros resultados. Para ello, usaremos la función trainControl() dentro del paquete caret. Esta se va a encargar de definir la configuración óptima del modelo. La función presenta tres parámetros:

  • method: “boot”, “boot632”, “optimism_boot”, “boot_all”, “cv”, “repeatedcv”, “LOOCV”,etc…

  • number: establece el numero de partes o bloques a dividir el conjunto de datos del mismo tamaño.

  • repeat: número de repeticiones.

La función selecciona el valor que da el mejor resultado.

fitControl <- trainControl(method = "repeatedcv",
                           number=5,
                           repeats=5)

Entrenamiento de un clasificador KNN (k-nearest neighbors)

La función train() del paquete caret es la que contiene la lógica para definir el clasificador de machine learning. Aspectos a destacar en el uso de la función:

  • clases ~ : Indica que se emplearán todos los atributos.

  • El parámetro data contiene las variables predictoras.

  • El parámetro method indica el métdo del clasificador a emplear, en este caso knn.

  • trControl establece como vamos a definir los parámetros del modelo, ver paso anterior.

En el caso de usar un clasificador KNN como método no paramétrico utilizado en clasificación, este predice o clasifica una muestra de datos utilizando la muestra más cercana a k de los datos de entrenamiento.

Por tanto, un clasificador KNN depende en gran medida de cómo se defina la distancia entre las muestras. Aunque hay muchas métricas de distancia, la distancia euclidiana es la que se utiliza habitualmente.

Dado que la distancia entre las muestras es crítica, se recomienda preprocesar (centrar y escalar) las variables predictoras antes de ejecutar el clasificador KNN. Esto elimina los sesgos y permite que todos los predictores sean tratados por igual al calcular la distancia.

Las ventajas del clasificador KNN son:

  • Es un clasificador sencillo que puede implementarse fácilmente.

  • Resulta apropiado para manejar clases multimodales.

Como inconveniente, requiere el cálculo de la distancia de los vecinos más cercanos, lo que puede demandar una alta capacidad de cómputo, sobre todo si el conjunto de datos de entrenamiento es grande.

knnFit <- train(leyenda ~ .,data=training,
               method="kknn",
               preProcess=c("center","scale"),
               trControl = fitControl)
print (knnFit)
## k-Nearest Neighbors 
## 
## 880 samples
##  10 predictor
##  13 classes: '--', 'Bosques Mixtos', 'Cadufolios', 'Castaños/Caducifolios', 'Eucalipto', 'Matorral', 'Pinos', 'Pinos/Pinsapos', 'Pinsapos', 'Quercineas', 'Ribera', 'Suelos?', 'Veg. Dispersa' 
## 
## Pre-processing: centered (10), scaled (10) 
## Resampling: Cross-Validated (5 fold, repeated 5 times) 
## Summary of sample sizes: 704, 704, 704, 704, 704, 704, ... 
## Resampling results across tuning parameters:
## 
##   kmax  Accuracy   Kappa    
##   5     0.3715909  0.2954465
##   7     0.3804545  0.3045002
##   9     0.3822727  0.3054572
## 
## Tuning parameter 'distance' was held constant at a value of 2
## Tuning
##  parameter 'kernel' was held constant at a value of optimal
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were kmax = 9, distance = 2 and kernel
##  = optimal.

Si graficamos los resultados podemos observar cual es el número adecuado de vecinos.

plot (knnFit)

Así como la configuaración del modelo

knnFit$finalModel
## 
## Call:
## kknn::train.kknn(formula = .outcome ~ ., data = dat, kmax = param$kmax,     distance = param$distance, kernel = as.character(param$kernel))
## 
## Type of response variable: nominal
## Minimal misclassification: 0.6079545
## Best kernel: optimal
## Best k: 7

Y la importancia de las variables empleadas a través de la función varImp().

knnvarImp <-varImp(knnFit,compete=FALSE)
plot(knnvarImp)

Además del entrenamiento, resulta necesario realizar un testeo del mismo, empleando para ello la función predict()

pred_knnFit <-predict(knnFit,newdata=testing)

Una vez realizado la fase de testeo sera necesario realizar un análisis de la calidad obtenida

confusionMatrix(data=pred_knnFit,testing$leyenda)
## Confusion Matrix and Statistics
## 
##                        Reference
## Prediction              -- Bosques Mixtos Cadufolios Castaños/Caducifolios
##   --                     3              2          0                     0
##   Bosques Mixtos         1             14          1                     0
##   Cadufolios             0              2          2                     1
##   Castaños/Caducifolios  0              0          0                     9
##   Eucalipto              1              2          0                     0
##   Matorral               1              3          1                     0
##   Pinos                  4              3          0                     0
##   Pinos/Pinsapos         0              0          0                     0
##   Pinsapos               0              1          0                     0
##   Quercineas             0              0          0                     0
##   Ribera                 0              1          0                     0
##   Suelos?                0              2          1                     0
##   Veg. Dispersa          0              0          5                     0
##                        Reference
## Prediction              Eucalipto Matorral Pinos Pinos/Pinsapos Pinsapos
##   --                            0        1     1              0        0
##   Bosques Mixtos                1        2     4              1        0
##   Cadufolios                    0        2     1              0        1
##   Castaños/Caducifolios         0        0     0              0        0
##   Eucalipto                     4        0     1              1        0
##   Matorral                      1        6     1              0        0
##   Pinos                         2        3    20              4        3
##   Pinos/Pinsapos                1        0     1              0        0
##   Pinsapos                      0        0     3              0        5
##   Quercineas                    0        0     1              4        1
##   Ribera                        1        0     2              0        0
##   Suelos?                       0        2     1              0        0
##   Veg. Dispersa                 0        4     4              0        0
##                        Reference
## Prediction              Quercineas Ribera Suelos? Veg. Dispersa
##   --                             0      1       1             0
##   Bosques Mixtos                 4      1       0             2
##   Cadufolios                     0      0       0             2
##   Castaños/Caducifolios          0      0       0             0
##   Eucalipto                      0      0       1             1
##   Matorral                       2      0       2             7
##   Pinos                          2      2       2             1
##   Pinos/Pinsapos                 1      1       0             1
##   Pinsapos                       0      1       0             0
##   Quercineas                     9      2       1             0
##   Ribera                         0      1       0             0
##   Suelos?                        0      0       0             2
##   Veg. Dispersa                  2      1       3            14
## 
## Overall Statistics
##                                           
##                Accuracy : 0.3955          
##                  95% CI : (0.3304, 0.4634)
##     No Information Rate : 0.1818          
##     P-Value [Acc > NIR] : 1.086e-13       
##                                           
##                   Kappa : 0.3214          
##                                           
##  Mcnemar's Test P-Value : NA              
## 
## Statistics by Class:
## 
##                      Class: -- Class: Bosques Mixtos Class: Cadufolios
## Sensitivity            0.30000               0.46667          0.200000
## Specificity            0.97143               0.91053          0.957143
## Pos Pred Value         0.33333               0.45161          0.181818
## Neg Pred Value         0.96682               0.91534          0.961722
## Prevalence             0.04545               0.13636          0.045455
## Detection Rate         0.01364               0.06364          0.009091
## Detection Prevalence   0.04091               0.14091          0.050000
## Balanced Accuracy      0.63571               0.68860          0.578571
##                      Class: Castaños/Caducifolios Class: Eucalipto
## Sensitivity                               0.90000          0.40000
## Specificity                               1.00000          0.96667
## Pos Pred Value                            1.00000          0.36364
## Neg Pred Value                            0.99526          0.97129
## Prevalence                                0.04545          0.04545
## Detection Rate                            0.04091          0.01818
## Detection Prevalence                      0.04091          0.05000
## Balanced Accuracy                         0.95000          0.68333
##                      Class: Matorral Class: Pinos Class: Pinos/Pinsapos
## Sensitivity                  0.30000      0.50000               0.00000
## Specificity                  0.91000      0.85556               0.97619
## Pos Pred Value               0.25000      0.43478               0.00000
## Neg Pred Value               0.92857      0.88506               0.95349
## Prevalence                   0.09091      0.18182               0.04545
## Detection Rate               0.02727      0.09091               0.00000
## Detection Prevalence         0.10909      0.20909               0.02273
## Balanced Accuracy            0.60500      0.67778               0.48810
##                      Class: Pinsapos Class: Quercineas Class: Ribera
## Sensitivity                  0.50000           0.45000      0.100000
## Specificity                  0.97619           0.95500      0.980952
## Pos Pred Value               0.50000           0.50000      0.200000
## Neg Pred Value               0.97619           0.94554      0.958140
## Prevalence                   0.04545           0.09091      0.045455
## Detection Rate               0.02273           0.04091      0.004545
## Detection Prevalence         0.04545           0.08182      0.022727
## Balanced Accuracy            0.73810           0.70250      0.540476
##                      Class: Suelos? Class: Veg. Dispersa
## Sensitivity                 0.00000              0.46667
## Specificity                 0.96190              0.90000
## Pos Pred Value              0.00000              0.42424
## Neg Pred Value              0.95283              0.91444
## Prevalence                  0.04545              0.13636
## Detection Rate              0.00000              0.06364
## Detection Prevalence        0.03636              0.15000
## Balanced Accuracy           0.48095              0.68333

Si los resultados en el control de calidad son positivos quedaría aplicar el modelo sobre la imagen para obtener la clasificación. (Este proceso requiere tiempo de computo dependiendo de la máquina empleada)

LC_knnFit <- predict(sentinel_o,knnFit)
mapview(LC_knnFit)
## Warning in rasterCheckSize(x, maxpixels = maxpixels): maximum number of pixels for Raster* viewing is 5e+05 ; 
## the supplied Raster* has 4220964 
##  ... decreasing Raster* resolution to 5e+05 pixels
##  to view full resolution set 'maxpixels =  4220964 '
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded ellps WGS 84 in Proj4 definition: +proj=merc +a=6378137
## +b=6378137 +lat_ts=0 +lon_0=0 +x_0=0 +y_0=0 +k=1 +units=m +nadgrids=@null
## +wktext +no_defs +type=crs
## Warning in showSRID(uprojargs, format = "PROJ", multiline = "NO", prefer_proj =
## prefer_proj): Discarded datum World Geodetic System 1984 in Proj4 definition

Entrenar a un clasificador ANN (Artificial Neural Networks)

Un clasificador ANN es un método basado en simular el funcionamiento del cerebro humano para: a) adquisición de conocimientos, b) recordar, c) sintetizar y, d) resolver problemas. Hay distintos clasificadore KNN com el perceptrón multicapa (MLP); mapas de características auto-organizados (SOM) de Kohonen; las redes de Hopfield; el clasificador de Carpenter/Grossberg.

De todos ellos, el de tipo MLP es uno de los modelos de red neuronal más utilizados, generalmente consta de tres o más capas:

  • Una capa de entrada: consta de uno o varios elementos de procesamiento que presentan los datos de entrenamiento.

  • Una o más capas ocultas, responsable de la representación interna de los datos, así como de la transformación de la información entre las capas de entrada y de salida.

  • Una capa de salida: consta de uno o varios elementos de procesamiento que almacenan los resultados de la red.

Las ventajas de un clasificadores KNN son que:

  • Los datos no deben estar sujetos a seguir una distribución normal.

  • Tienen la capacidad de generar límites de decisión no lineales.

  • Capacidad de aprender patrones complejos.

Sin embargo, los clasificadores de KNN son propensos al overfiting sobreajuste, siendo complejo el diseñar una red neuronal eficaz. El overfitting se debe a que el modelo se ajustará a aprender los casos particulares que le enseñamos, siendo incapaz de reconocer nuevos datos de entrada.

annFit <- train(leyenda ~ ., data = training,
                method = "nnet",
                preProcess = c("center", "scale"),
                trControl = fitControl)
print (annFit)
## Neural Network 
## 
## 880 samples
##  10 predictor
##  13 classes: '--', 'Bosques Mixtos', 'Cadufolios', 'Castaños/Caducifolios', 'Eucalipto', 'Matorral', 'Pinos', 'Pinos/Pinsapos', 'Pinsapos', 'Quercineas', 'Ribera', 'Suelos?', 'Veg. Dispersa' 
## 
## Pre-processing: centered (10), scaled (10) 
## Resampling: Cross-Validated (5 fold, repeated 5 times) 
## Summary of sample sizes: 704, 704, 704, 704, 704, 704, ... 
## Resampling results across tuning parameters:
## 
##   size  decay  Accuracy   Kappa    
##   1     0e+00  0.2854545  0.1586734
##   1     1e-04  0.2815909  0.1543827
##   1     1e-01  0.2643182  0.1236267
##   3     0e+00  0.3820455  0.2946602
##   3     1e-04  0.3943182  0.3087801
##   3     1e-01  0.3850000  0.2908279
##   5     0e+00  0.4190909  0.3427078
##   5     1e-04  0.4245455  0.3486156
##   5     1e-01  0.4234091  0.3399757
## 
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were size = 5 and decay = 1e-04.

La siguiente figura muestra la relación entre pesos y capas ocultas de la red. En este ejemplo la mejor configuración del modelo ANN se obtiene con 5 nodos y un peso igual a 0.1.

plot (annFit)

Así, la configuración del modelo ANN es: - Se ha entrenado con 10 variables (las 10 bandas espectrales de Sentinel 2) - 5 nodos en la capa oculta de la red - 9 clase para la capa de salida

annFit$finalModel
## a 10-5-13 network with 133 weights
## inputs: sentinel_o_b01 sentinel_o_b02 sentinel_o_b03 sentinel_o_b04 sentinel_o_b05 sentinel_o_b06 sentinel_o_b07 sentinel_o_b08 sentinel_o_b09 sentinel_o_b10 
## output(s): .outcome 
## options were - softmax modelling  decay=1e-04

Si deseamos visualizar la red emplearemos la función plotnet.

plotnet(annFit$finalModel)

Y mediante la función olden podemos analizar a importancia relativa de las variables predictoras. En este caso la banda 1 es la que más importancia presenta, todo lo contrario que la banda 5.

olden(annFit)

Y ahora, realizamos la predicción y el control de calidad sobre esta.

pred_annFit<- predict(annFit, newdata = testing)
confusionMatrix(data = pred_annFit, testing$leyenda)
## Confusion Matrix and Statistics
## 
##                        Reference
## Prediction              -- Bosques Mixtos Cadufolios Castaños/Caducifolios
##   --                     2              0          0                     0
##   Bosques Mixtos         2             15          0                     0
##   Cadufolios             0              1          1                     0
##   Castaños/Caducifolios  0              0          0                     9
##   Eucalipto              0              1          0                     0
##   Matorral               0              1          2                     0
##   Pinos                  4              7          1                     0
##   Pinos/Pinsapos         0              0          0                     0
##   Pinsapos               0              0          0                     0
##   Quercineas             0              1          1                     0
##   Ribera                 0              0          1                     0
##   Suelos?                0              0          0                     0
##   Veg. Dispersa          2              4          4                     1
##                        Reference
## Prediction              Eucalipto Matorral Pinos Pinos/Pinsapos Pinsapos
##   --                            0        3     1              0        0
##   Bosques Mixtos                4        6     6              3        0
##   Cadufolios                    0        0     1              0        1
##   Castaños/Caducifolios         0        0     0              0        0
##   Eucalipto                     5        0     1              0        0
##   Matorral                      0        6     0              0        0
##   Pinos                         0        1    28              6        5
##   Pinos/Pinsapos                0        0     0              0        0
##   Pinsapos                      0        2     0              0        2
##   Quercineas                    0        0     0              0        0
##   Ribera                        1        0     0              0        0
##   Suelos?                       0        0     0              0        0
##   Veg. Dispersa                 0        2     3              1        2
##                        Reference
## Prediction              Quercineas Ribera Suelos? Veg. Dispersa
##   --                             1      0       0             2
##   Bosques Mixtos                 3      4       1             2
##   Cadufolios                     0      0       0             0
##   Castaños/Caducifolios          0      0       0             1
##   Eucalipto                      1      1       0             0
##   Matorral                       0      0       2             8
##   Pinos                          1      3       5             2
##   Pinos/Pinsapos                 0      0       0             0
##   Pinsapos                       0      0       0             0
##   Quercineas                    10      1       0             3
##   Ribera                         0      0       0             0
##   Suelos?                        0      0       0             2
##   Veg. Dispersa                  4      1       2            10
## 
## Overall Statistics
##                                          
##                Accuracy : 0.4            
##                  95% CI : (0.3347, 0.468)
##     No Information Rate : 0.1818         
##     P-Value [Acc > NIR] : 3.613e-14      
##                                          
##                   Kappa : 0.3138         
##                                          
##  Mcnemar's Test P-Value : NA             
## 
## Statistics by Class:
## 
##                      Class: -- Class: Bosques Mixtos Class: Cadufolios
## Sensitivity           0.200000               0.50000          0.100000
## Specificity           0.966667               0.83684          0.985714
## Pos Pred Value        0.222222               0.32609          0.250000
## Neg Pred Value        0.962085               0.91379          0.958333
## Prevalence            0.045455               0.13636          0.045455
## Detection Rate        0.009091               0.06818          0.004545
## Detection Prevalence  0.040909               0.20909          0.018182
## Balanced Accuracy     0.583333               0.66842          0.542857
##                      Class: Castaños/Caducifolios Class: Eucalipto
## Sensitivity                               0.90000          0.50000
## Specificity                               0.99524          0.98095
## Pos Pred Value                            0.90000          0.55556
## Neg Pred Value                            0.99524          0.97630
## Prevalence                                0.04545          0.04545
## Detection Rate                            0.04091          0.02273
## Detection Prevalence                      0.04545          0.04091
## Balanced Accuracy                         0.94762          0.74048
##                      Class: Matorral Class: Pinos Class: Pinos/Pinsapos
## Sensitivity                  0.30000       0.7000               0.00000
## Specificity                  0.93500       0.8056               1.00000
## Pos Pred Value               0.31579       0.4444                   NaN
## Neg Pred Value               0.93035       0.9236               0.95455
## Prevalence                   0.09091       0.1818               0.04545
## Detection Rate               0.02727       0.1273               0.00000
## Detection Prevalence         0.08636       0.2864               0.00000
## Balanced Accuracy            0.61750       0.7528               0.50000
##                      Class: Pinsapos Class: Quercineas Class: Ribera
## Sensitivity                 0.200000           0.50000      0.000000
## Specificity                 0.990476           0.97000      0.990476
## Pos Pred Value              0.500000           0.62500      0.000000
## Neg Pred Value              0.962963           0.95098      0.954128
## Prevalence                  0.045455           0.09091      0.045455
## Detection Rate              0.009091           0.04545      0.000000
## Detection Prevalence        0.018182           0.07273      0.009091
## Balanced Accuracy           0.595238           0.73500      0.495238
##                      Class: Suelos? Class: Veg. Dispersa
## Sensitivity                0.000000              0.33333
## Specificity                0.990476              0.86316
## Pos Pred Value             0.000000              0.27778
## Neg Pred Value             0.954128              0.89130
## Prevalence                 0.045455              0.13636
## Detection Rate             0.000000              0.04545
## Detection Prevalence       0.009091              0.16364
## Balanced Accuracy          0.495238              0.59825

Entrenar a un clasificador SVM (Support Vector Machine)

Un clasificador SVM se basa en modelos estadísticos los cuales, con caracter tiene por objeto localizar una frontera de decisión (separación) óptima que maximice el margen entre dos clases.

La ubicación del límite de decisión dependerá solo de un subconjunto de puntos de datos de entrenamiento que están más cerca de él. Este subconjunto de puntos de datos de entrenamiento más cercanos al límite de decisión se conoce como vectores de soporte. A lo largo del tiempo han aparecido evolutivos incorporando funciones de coste, funciones para permitir límites de clase no lineales.

Además, funciones polinómicas, base radial o tangente hiperbólica, se desarrollaron para transformar el conjunto de datos de entrenamiento en un espacio de características de mayor dimensión para los problemas de clasificación no lineal.

Las ventajas de un clasificador SVM son:

  • Se elimina el problema de la posible suposición de distribución normal de los datos.

  • Uso de una funciones para resolver problemas complejos.

  • Se adapta relativamente bien a datos de alta dimensión.

No obstante, como inconventiente, este clasificador necesita de más tiempo de entrenamiento.

svm_model<-train(leyenda~.,data=training,
                  method = "svmRadial",
                  trControl = fitControl,
                  preProc = c("center", "scale"),
                  tuneLength = 3)

Al hacer una impresión del modelo se observa como este tiene una calidad igual a 0.36 y un coste igual a 1.

print (svm_model)
## Support Vector Machines with Radial Basis Function Kernel 
## 
## 880 samples
##  10 predictor
##  13 classes: '--', 'Bosques Mixtos', 'Cadufolios', 'Castaños/Caducifolios', 'Eucalipto', 'Matorral', 'Pinos', 'Pinos/Pinsapos', 'Pinsapos', 'Quercineas', 'Ribera', 'Suelos?', 'Veg. Dispersa' 
## 
## Pre-processing: centered (10), scaled (10) 
## Resampling: Cross-Validated (5 fold, repeated 5 times) 
## Summary of sample sizes: 704, 704, 704, 704, 704, 704, ... 
## Resampling results across tuning parameters:
## 
##   C     Accuracy   Kappa    
##   0.25  0.3131818  0.1843286
##   0.50  0.3295455  0.2112172
##   1.00  0.3665909  0.2671307
## 
## Tuning parameter 'sigma' was held constant at a value of 0.5696137
## Accuracy was used to select the optimal model using the largest value.
## The final values used for the model were sigma = 0.5696137 and C = 1.
plot (svm_model)

El siguiente paso es mostrar la importancia de las variables de entrada

svm_varImp <- varImp(svm_model, compete = FALSE)
plot(svm_varImp)

Y finalmente realizar un control de calidad

pred_svm<- predict(svm_model, newdata = testing)
confusionMatrix(data = pred_svm, testing$leyenda)
## Confusion Matrix and Statistics
## 
##                        Reference
## Prediction              -- Bosques Mixtos Cadufolios Castaños/Caducifolios
##   --                     3              0          0                     0
##   Bosques Mixtos         1              5          0                     0
##   Cadufolios             0              2          0                     2
##   Castaños/Caducifolios  0              0          0                     8
##   Eucalipto              0              0          1                     0
##   Matorral               0              3          2                     0
##   Pinos                  6             16          1                     0
##   Pinos/Pinsapos         0              0          0                     0
##   Pinsapos               0              1          1                     0
##   Quercineas             0              2          0                     0
##   Ribera                 0              0          0                     0
##   Suelos?                0              0          0                     0
##   Veg. Dispersa          0              1          5                     0
##                        Reference
## Prediction              Eucalipto Matorral Pinos Pinos/Pinsapos Pinsapos
##   --                            0        2     0              0        0
##   Bosques Mixtos                1        4     2              2        0
##   Cadufolios                    0        1     0              0        0
##   Castaños/Caducifolios         0        0     0              0        0
##   Eucalipto                     0        0     1              0        0
##   Matorral                      0       10     0              0        1
##   Pinos                         6        2    33              8        6
##   Pinos/Pinsapos                0        0     0              0        0
##   Pinsapos                      2        0     1              0        2
##   Quercineas                    0        0     0              0        1
##   Ribera                        0        0     0              0        0
##   Suelos?                       0        0     0              0        0
##   Veg. Dispersa                 1        1     3              0        0
##                        Reference
## Prediction              Quercineas Ribera Suelos? Veg. Dispersa
##   --                             0      0       0             0
##   Bosques Mixtos                 2      0       1             3
##   Cadufolios                     0      0       0             0
##   Castaños/Caducifolios          0      0       0             0
##   Eucalipto                      0      0       1             0
##   Matorral                       0      0       1             5
##   Pinos                          8      7       3             2
##   Pinos/Pinsapos                 0      0       0             0
##   Pinsapos                       0      0       0             0
##   Quercineas                     7      3       0             2
##   Ribera                         0      0       0             0
##   Suelos?                        0      0       0             0
##   Veg. Dispersa                  3      0       4            18
## 
## Overall Statistics
##                                          
##                Accuracy : 0.3909         
##                  95% CI : (0.326, 0.4588)
##     No Information Rate : 0.1818         
##     P-Value [Acc > NIR] : 3.202e-13      
##                                          
##                   Kappa : 0.2939         
##                                          
##  Mcnemar's Test P-Value : NA             
## 
## Statistics by Class:
## 
##                      Class: -- Class: Bosques Mixtos Class: Cadufolios
## Sensitivity            0.30000               0.16667           0.00000
## Specificity            0.99048               0.91579           0.97619
## Pos Pred Value         0.60000               0.23810           0.00000
## Neg Pred Value         0.96744               0.87437           0.95349
## Prevalence             0.04545               0.13636           0.04545
## Detection Rate         0.01364               0.02273           0.00000
## Detection Prevalence   0.02273               0.09545           0.02273
## Balanced Accuracy      0.64524               0.54123           0.48810
##                      Class: Castaños/Caducifolios Class: Eucalipto
## Sensitivity                               0.80000          0.00000
## Specificity                               1.00000          0.98571
## Pos Pred Value                            1.00000          0.00000
## Neg Pred Value                            0.99057          0.95392
## Prevalence                                0.04545          0.04545
## Detection Rate                            0.03636          0.00000
## Detection Prevalence                      0.03636          0.01364
## Balanced Accuracy                         0.90000          0.49286
##                      Class: Matorral Class: Pinos Class: Pinos/Pinsapos
## Sensitivity                  0.50000       0.8250               0.00000
## Specificity                  0.94000       0.6389               1.00000
## Pos Pred Value               0.45455       0.3367                   NaN
## Neg Pred Value               0.94949       0.9426               0.95455
## Prevalence                   0.09091       0.1818               0.04545
## Detection Rate               0.04545       0.1500               0.00000
## Detection Prevalence         0.10000       0.4455               0.00000
## Balanced Accuracy            0.72000       0.7319               0.50000
##                      Class: Pinsapos Class: Quercineas Class: Ribera
## Sensitivity                 0.200000           0.35000       0.00000
## Specificity                 0.976190           0.96000       1.00000
## Pos Pred Value              0.285714           0.46667           NaN
## Neg Pred Value              0.962441           0.93659       0.95455
## Prevalence                  0.045455           0.09091       0.04545
## Detection Rate              0.009091           0.03182       0.00000
## Detection Prevalence        0.031818           0.06818       0.00000
## Balanced Accuracy           0.588095           0.65500       0.50000
##                      Class: Suelos? Class: Veg. Dispersa
## Sensitivity                 0.00000              0.60000
## Specificity                 1.00000              0.90526
## Pos Pred Value                  NaN              0.50000
## Neg Pred Value              0.95455              0.93478
## Prevalence                  0.04545              0.13636
## Detection Rate              0.00000              0.08182
## Detection Prevalence        0.00000              0.16364
## Balanced Accuracy           0.50000              0.75263

Entrenar a un clasificador RF (Random Forest)

El clasificador RF es un método de aprendizaje automático de conjunto, que utiliza el muestreo bootstrap para construir muchos modelos de árboles de decisión individuales.

Usa un subconjunto aleatorio de variables predictoras (por ejemplo, las bandas Sentinel) para dividir los datos de observación en subconjuntos homogéneos, que se utilizan para construir cada modelo de árbol de decisión y una predicción. Luego, se promedian las predicciones del modelo de árbol de decisión individual para producir el etiquetado final.

El clasificador RF se ha utilizado con éxito para la clasificación de imágenes de teledetección porque presentan estas ventajas:

  • Permiten manejar grandes cantidades de datos.

  • Están libres de supuestos de distribución normal.

  • Son robustos a los valores atípicos y al ruido

Sin embargo, como inconventientes:

  • No es fácil interpretar los resultados del modelo RF.

  • Esta sesgado a favor de las variables predictoras con muchos niveles de categorías diferentes.

rf_model<-train(leyenda~.,data=training, method="rf",
                trControl=fitControl,
                 prox=TRUE,
                 fitBest = FALSE,
                 returnData = TRUE)
print(rf_model)
## Random Forest 
## 
## 880 samples
##  10 predictor
##  13 classes: '--', 'Bosques Mixtos', 'Cadufolios', 'Castaños/Caducifolios', 'Eucalipto', 'Matorral', 'Pinos', 'Pinos/Pinsapos', 'Pinsapos', 'Quercineas', 'Ribera', 'Suelos?', 'Veg. Dispersa' 
## 
## No pre-processing
## Resampling: Cross-Validated (5 fold, repeated 5 times) 
## Summary of sample sizes: 704, 704, 704, 704, 704, 704, ... 
## Resampling results across tuning parameters:
## 
##   mtry  Accuracy   Kappa    
##    2    0.4195455  0.3434671
##    6    0.4284091  0.3542228
##   10    0.4304545  0.3566526
## 
## Accuracy was used to select the optimal model using the largest value.
## The final value used for the model was mtry = 10.
plot(rf_model)

rf_model$finalModel
## 
## Call:
##  randomForest(x = x, y = y, mtry = min(param$mtry, ncol(x)), proximity = TRUE,      fitBest = FALSE, returnData = TRUE) 
##                Type of random forest: classification
##                      Number of trees: 500
## No. of variables tried at each split: 10
## 
##         OOB estimate of  error rate: 56.82%
## Confusion matrix:
##                       -- Bosques Mixtos Cadufolios Castaños/Caducifolios
## --                    24              4          0                     0
## Bosques Mixtos         7             51          1                     0
## Cadufolios             0              1         16                     1
## Castaños/Caducifolios  0              1          1                    36
## Eucalipto              0              6          0                     0
## Matorral               3              7          4                     1
## Pinos                  1             25          3                     0
## Pinos/Pinsapos         1              6          1                     0
## Pinsapos               0              2          0                     0
## Quercineas             2             10          1                     0
## Ribera                 0              5          0                     0
## Suelos?                0              5          4                     0
## Veg. Dispersa          0              9         12                     4
##                       Eucalipto Matorral Pinos Pinos/Pinsapos Pinsapos
## --                            0        1     7              3        0
## Bosques Mixtos                1        4    26              0        3
## Cadufolios                    2        2     3              0        0
## Castaños/Caducifolios         0        0     0              0        0
## Eucalipto                    11        1    12              0        0
## Matorral                      0       26     9              1        0
## Pinos                         3        4    93              6        6
## Pinos/Pinsapos                0        1    15              5        4
## Pinsapos                      1        1    11              1       20
## Quercineas                    1        0    18              0        1
## Ribera                        2        3     6              0        0
## Suelos?                       3        6     4              1        3
## Veg. Dispersa                 5       18    14              1        0
##                       Quercineas Ribera Suelos? Veg. Dispersa class.error
## --                             0      0       0             1   0.4000000
## Bosques Mixtos                12      0       2            13   0.5750000
## Cadufolios                     0      0       0            15   0.6000000
## Castaños/Caducifolios          0      0       0             2   0.1000000
## Eucalipto                      0      4       2             4   0.7250000
## Matorral                       2      1       2            24   0.6750000
## Pinos                          9      2       1             7   0.4187500
## Pinos/Pinsapos                 3      0       0             4   0.8750000
## Pinsapos                       0      1       2             1   0.5000000
## Quercineas                    34      3       1             9   0.5750000
## Ribera                         4     16       0             4   0.6000000
## Suelos?                        0      2       1            11   0.9750000
## Veg. Dispersa                  6      1       3            47   0.6083333
rf_varImp <- varImp(rf_model, compete = FALSE)
plot(rf_varImp)

Como en los clasificadores estudiados anteriormente realizaremos un control de calidad.

pred_rf <- predict(rf_model$finalModel,
            newdata = testing)
confusionMatrix(data = pred_rf, testing$leyenda)
## Confusion Matrix and Statistics
## 
##                        Reference
## Prediction              -- Bosques Mixtos Cadufolios Castaños/Caducifolios
##   --                     3              2          0                     0
##   Bosques Mixtos         1             13          0                     0
##   Cadufolios             0              2          3                     1
##   Castaños/Caducifolios  0              0          0                     9
##   Eucalipto              0              1          1                     0
##   Matorral               0              2          1                     0
##   Pinos                  6              6          0                     0
##   Pinos/Pinsapos         0              0          0                     0
##   Pinsapos               0              1          0                     0
##   Quercineas             0              2          0                     0
##   Ribera                 0              0          1                     0
##   Suelos?                0              0          0                     0
##   Veg. Dispersa          0              1          4                     0
##                        Reference
## Prediction              Eucalipto Matorral Pinos Pinos/Pinsapos Pinsapos
##   --                            2        2     1              0        0
##   Bosques Mixtos                1        2     1              1        0
##   Cadufolios                    0        2     1              0        0
##   Castaños/Caducifolios         0        0     0              0        0
##   Eucalipto                     1        0     2              1        1
##   Matorral                      0        8     0              0        0
##   Pinos                         3        1    22              6        5
##   Pinos/Pinsapos                0        1     1              0        0
##   Pinsapos                      0        0     2              1        2
##   Quercineas                    0        0     3              0        0
##   Ribera                        1        0     0              0        0
##   Suelos?                       0        0     0              0        1
##   Veg. Dispersa                 2        4     7              1        1
##                        Reference
## Prediction              Quercineas Ribera Suelos? Veg. Dispersa
##   --                             0      0       1             0
##   Bosques Mixtos                 6      3       0             1
##   Cadufolios                     0      0       1             2
##   Castaños/Caducifolios          0      0       0             0
##   Eucalipto                      0      0       0             0
##   Matorral                       1      0       2             5
##   Pinos                          4      4       3             2
##   Pinos/Pinsapos                 0      0       0             1
##   Pinsapos                       0      0       0             0
##   Quercineas                     7      1       0             0
##   Ribera                         0      2       0             0
##   Suelos?                        0      0       1             3
##   Veg. Dispersa                  2      0       2            16
## 
## Overall Statistics
##                                           
##                Accuracy : 0.3955          
##                  95% CI : (0.3304, 0.4634)
##     No Information Rate : 0.1818          
##     P-Value [Acc > NIR] : 1.086e-13       
##                                           
##                   Kappa : 0.3138          
##                                           
##  Mcnemar's Test P-Value : NA              
## 
## Statistics by Class:
## 
##                      Class: -- Class: Bosques Mixtos Class: Cadufolios
## Sensitivity            0.30000               0.43333           0.30000
## Specificity            0.96190               0.91579           0.95714
## Pos Pred Value         0.27273               0.44828           0.25000
## Neg Pred Value         0.96651               0.91099           0.96635
## Prevalence             0.04545               0.13636           0.04545
## Detection Rate         0.01364               0.05909           0.01364
## Detection Prevalence   0.05000               0.13182           0.05455
## Balanced Accuracy      0.63095               0.67456           0.62857
##                      Class: Castaños/Caducifolios Class: Eucalipto
## Sensitivity                               0.90000         0.100000
## Specificity                               1.00000         0.971429
## Pos Pred Value                            1.00000         0.142857
## Neg Pred Value                            0.99526         0.957746
## Prevalence                                0.04545         0.045455
## Detection Rate                            0.04091         0.004545
## Detection Prevalence                      0.04091         0.031818
## Balanced Accuracy                         0.95000         0.535714
##                      Class: Matorral Class: Pinos Class: Pinos/Pinsapos
## Sensitivity                  0.40000       0.5500               0.00000
## Specificity                  0.94500       0.7778               0.98571
## Pos Pred Value               0.42105       0.3548               0.00000
## Neg Pred Value               0.94030       0.8861               0.95392
## Prevalence                   0.09091       0.1818               0.04545
## Detection Rate               0.03636       0.1000               0.00000
## Detection Prevalence         0.08636       0.2818               0.01364
## Balanced Accuracy            0.67250       0.6639               0.49286
##                      Class: Pinsapos Class: Quercineas Class: Ribera
## Sensitivity                 0.200000           0.35000      0.200000
## Specificity                 0.980952           0.97000      0.990476
## Pos Pred Value              0.333333           0.53846      0.500000
## Neg Pred Value              0.962617           0.93720      0.962963
## Prevalence                  0.045455           0.09091      0.045455
## Detection Rate              0.009091           0.03182      0.009091
## Detection Prevalence        0.027273           0.05909      0.018182
## Balanced Accuracy           0.590476           0.66000      0.595238
##                      Class: Suelos? Class: Veg. Dispersa
## Sensitivity                0.100000              0.53333
## Specificity                0.980952              0.87368
## Pos Pred Value             0.200000              0.40000
## Neg Pred Value             0.958140              0.92222
## Prevalence                 0.045455              0.13636
## Detection Rate             0.004545              0.07273
## Detection Prevalence       0.022727              0.18182
## Balanced Accuracy          0.540476              0.70351

Comparacón de los clasificadores estudiados

La comparación se realizará siguiendo una correlación cruzada, empleando para esto la función resample().

resamps <- resamples(list(knn = knnFit,
                          ann = annFit,
                          svm = svm_model,
                          rf = rf_model))

Representamos mediante un boxplot los valores de accuracy y kappa para realizar una evaluación de la calidad entre clasificadores empleados.

bwplot(resamps, layout = c(3, 1))

A continuación se realizará la predicción de cada uno de los clasificadores anteriormente desarrollados.

LC_knnFit <-predict(sentinel_o,knnFit)
LC_ann <-predict(sentinel_o,annFit)
LC_svm <-predict(sentinel_o,svm_model)
LC_rf <-predict(sentinel_o,rf_model)

Clasificación digital empleando datos multitemporales

El objetivo es obtener una clasificación de la misma zona de estudio que en el tutorial 1 pero en este caso considerando datos multitemporales correspondientes a tres escenas de Sentinel 2: primeravera, verano y otoño.

Vamos en primer lugar a generar un rasterstack según fecha.

dir_in='./Material_practicas/Sentinel/'
archivos <-list.files(paste(dir_in,"o",sep="")
                      ,pattern = ".tif",
                      full.names = TRUE)
sentinel_o <- stack(archivos)

archivos <-list.files(paste(dir_in,"p",sep="")
                      ,pattern = ".tif",
                      full.names = TRUE)
sentinel_p <- stack(archivos)

archivos <-list.files(paste(dir_in,"v",sep="")
                      ,pattern = ".tif",
                      full.names = TRUE)
sentinel_v <- stack(archivos)
plotRGB(sentinel_o,r=7,g=2,b=3,stretch="lin")

plotRGB(sentinel_p,r=7,g=2,b=3,stretch="lin")

plotRGB(sentinel_v,r=7,g=2,b=3,stretch="lin")

Crearemos a continuación un rasterstack resultante de las escenas Sentinel 2 de las tres fechas, teniendo un total de 30 bandas espectrales (10 por cada una de las fechas)

sentinel <- stack(sentinel_o,sentinel_p,sentinel_v)

A partir de aquí, el desarrollo del clasificador será exactamente igual que en el caso de trabajar con una imagen, salvo que el tiempo de procesado será mayor